# below was the procedure to initially process the files; I have saved the resulting file, negating the need to re-run

datapath = "/Users/bjpessman/Documents/phd_research_code/Vibratory_Noise/recording_files"
setwd(datapath)
txt_files_ls = list.files(path = datapath, pattern = "*.txt")
txt_files_df <- lapply(txt_files_ls, function(x) 
    {read.table(file = x, header = TRUE, sep ="\t")})
combined_dfl <- do.call("rbind", lapply(txt_files_df, as.data.frame))
combined_dfl$Time <- strptime(combined_dfl$Time, format = "%H:%M:%S")
combined_dfl$Hour <-hour(combined_dfl$Time)
pca_results <- readRDS("/Users/bjpessman/Documents/phd_research_code/Vibratory_Noise/data/pca_results.rds") %>% 
  mutate(Site = site, .keep = "unused")
combined_dfl <- full_join(combined_dfl, pca_results, by = "Site")

# Calculate the daily average leq for each 24-hour recording at each site
dayavgl <- combined_dfl %>% 
  mutate(Substrate = fct_recode(factor(Substrate), "Plant" = "Plant", "Manmade" = "Artificial", "Manmade" = "Artifical")) %>% 
  rename(Traffic_Impact = `Traffic Impact`) %>% 
  group_by(Site, Visit, Mic, Category, Substrate, Material, Dim.1, Traffic_Impact) %>% 
  summarize(mean_leq = mean(Leq, na.rm = TRUE), 
            mean_aggent = mean(AggEnt, na.rm = TRUE))

prephouravgcatl <- combined_dfl %>% 
  group_by(Site, Visit, Mic, Hour, Category, Substrate, Material) %>% 
  summarize(meanleq = mean(Leq, na.rm = TRUE))

houravgcatl <- prephouravgcatl %>% 
  mutate(dark_light = as.numeric(ifelse(Hour >= 7, Hour - 7, Hour + 17)))

# get the average leq for each site
siteavgl <- combined_dfl %>% 
  group_by(Site, Category) %>% 
  summarize(mean_leq = mean(Leq)) 

saveRDS(dayavgl, "/Users/bjpessman/Documents/phd_research_code/Vibratory_Noise/data/dayavgl.rds")
saveRDS(houravgcatl, "/Users/bjpessman/Documents/phd_research_code/Vibratory_Noise/data/houravgcatl.rds")
saveRDS(siteavgl, "/Users/bjpessman/Documents/phd_research_code/Vibratory_Noise/data/siteavgl.rds")

Fun with Maps

## [1] "Urban" "Rural"
## quartz_off_screen 
##                 2

Sites are randomly distributed across Lincoln, Nebraska (urban) and into the surrounding rural area (rural).

## quartz_off_screen 
##                 2

Site averages vary from -55 dB to -70 dB (~15 dB difference). The sites with the highest Leq appear to occur near highly traveled roads (i.e., highways/interstates).

Spatial Variation in Vibratory Noise

Does vibratory noise vary over space?

## [1] 179 276

## [1] 179 276

## [1] 179 276

## [1] 179 276
## 
## Call:
## lm(formula = mean_leq ~ Dim.1 * Substrate, data = dayavgl)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7407 -1.8525 -0.4593  1.2537 10.5671 
## 
## Coefficients:
##                        Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)          -62.946760   0.257397 -244.552  < 2e-16 ***
## Dim.1                  1.738579   0.150514   11.551  < 2e-16 ***
## SubstratePlant        -1.285796   0.346281   -3.713 0.000245 ***
## Dim.1:SubstratePlant   0.004224   0.203036    0.021 0.983417    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.957 on 291 degrees of freedom
## Multiple R-squared:  0.5171, Adjusted R-squared:  0.5121 
## F-statistic: 103.9 on 3 and 291 DF,  p-value: < 2.2e-16

## quartz_off_screen 
##                 2

Daily average Leq has a significant positive relationship with Principal Component 1 - road vibratory noise (t = 10.257, df = 291, P < 0.001, adj.R^2 = 0.45). Daily average Leq was significantly higher on manmade material than plant material (t = -3.298, df = 291, P = 0.001). There is no interaction (t = -0.064, p = 0.949)

Let’s take a closer look at the substrate. Manmade structures - Paneling, Metal, Concrete, Brick, Wood Plant structures - Herb, Tree, Shrub, Vine

Brick, paneling, and shrubs carried the highest amplitude vibrations Bricks and herbs have the steepest slopes, which might suggest these substrates are affected by vibratory noise. Wood in quiet areas have high vibrations, probably as a result of people and pets walking on porches.

Temporal Variation in Vibratory Noise

Season

## Linear mixed model fit by REML ['lmerMod']
## Formula: mean_leq ~ Visit * Category + (1 | Site)
##    Data: dayavgl_20
## 
## REML criterion at convergence: 1178
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3807 -0.6166 -0.0366  0.4206  3.8739 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 5.743    2.397   
##  Residual             3.943    1.986   
## Number of obs: 268, groups:  Site, 21
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)          -62.4301     0.7051 -88.543
## Visit2                -0.4990     0.4490  -1.111
## Visit3                 0.1681     0.4371   0.385
## Visit4                -0.2262     0.4299  -0.526
## CategoryRural         -6.1084     1.3005  -4.697
## Visit2:CategoryRural   0.4400     0.7848   0.561
## Visit3:CategoryRural   1.3976     0.8142   1.717
## Visit4:CategoryRural   0.9722     0.7748   1.255
## 
## Correlation of Fixed Effects:
##             (Intr) Visit2 Visit3 Visit4 CtgryR Vs2:CR Vs3:CR
## Visit2      -0.360                                          
## Visit3      -0.366  0.574                                   
## Visit4      -0.379  0.593  0.604                            
## CategoryRrl -0.542  0.195  0.199  0.205                     
## Vst2:CtgryR  0.206 -0.572 -0.329 -0.339 -0.343              
## Vst3:CtgryR  0.197 -0.308 -0.537 -0.324 -0.330  0.545       
## Vst4:CtgryR  0.210 -0.329 -0.335 -0.555 -0.348  0.575  0.553
## $`emmeans of Visit`
##  Visit emmean    SE   df lower.CL upper.CL
##  1      -65.5 0.650 27.1    -66.8    -64.2
##  2      -65.8 0.634 24.5    -67.1    -64.5
##  3      -64.6 0.643 25.9    -65.9    -63.3
##  4      -65.2 0.631 24.0    -66.5    -63.9
## 
## Results are averaged over the levels of: Category 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $`pairwise differences of Visit`
##  1               estimate    SE  df t.ratio p.value
##  Visit1 - Visit2    0.279 0.392 242   0.711  0.8927
##  Visit1 - Visit3   -0.867 0.407 242  -2.129  0.1468
##  Visit1 - Visit4   -0.260 0.387 242  -0.671  0.9081
##  Visit2 - Visit3   -1.146 0.382 241  -3.002  0.0156
##  Visit2 - Visit4   -0.539 0.360 241  -1.498  0.4401
##  Visit3 - Visit4    0.607 0.376 241   1.614  0.3728
## 
## Results are averaged over the levels of: Category 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 4 estimates
##            R2m       R2c
## [1,] 0.3928398 0.7528386

## quartz_off_screen 
##                 2

To investigate whether noise varied across the season, we used a linear mixed model with visit number and category, and their interaction with site as a random factor. There was a trend that daily average Leq varied across the 2020 season (Chisq = 7.362, df = 265, P = 0.061, Conditional R^2 = 0.753). A post hoc test suggests that visit 3 was significantly louder than visit 2 (t = -3.002, P = 0.016). Also, urban areas are louder than rural areas (Chisq = 20.494, P < 0.001). There is no interaction between visit and category (Chisq = 3.493, P = 0.322).

Let’s look at date rather than visit.

## 
## Call:
## lm(formula = mean_leq ~ Day + Category, data = dayavgl_20)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1056 -2.0208 -0.6313  1.3565 12.0737 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -64.996920   2.124707 -30.591   <2e-16 ***
## Day             0.009532   0.008037   1.186    0.237    
## CategoryRural  -5.542657   0.413734 -13.397   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.056 on 265 degrees of freedom
## Multiple R-squared:  0.404,  Adjusted R-squared:  0.3995 
## F-statistic:  89.8 on 2 and 265 DF,  p-value: < 2.2e-16
## quartz_off_screen 
##                 2

We see similar results, but with no difference over time (F = 1.404, df = 264, P = 0.237, Adj. R^2 = 0.398).

## Linear mixed model fit by REML ['lmerMod']
## Formula: mean_leq ~ mean_harvest + (1 | Site)
##    Data: dayavgl_20_rural
## 
## REML criterion at convergence: 344.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7332 -0.5890 -0.1364  0.2263  3.9614 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 2.247    1.499   
##  Residual             4.038    2.010   
## Number of obs: 78, groups:  Site, 6
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)  -68.70767    0.73442 -93.554
## mean_harvest   0.07204    0.03397   2.121
## 
## Correlation of Fixed Effects:
##             (Intr)
## mean_harvst -0.456
##             R2m       R2c
## [1,] 0.03639155 0.3808811

## quartz_off_screen 
##                 2

We used USDA data on week end percent harvest in 2020 for field crops in Nebraska. This gave details on oats, wheat, dry beans, sorghum, corn, and soybeans. We restricted this list to corn and soybeans, as these are the major crops grown and harvested in Lancaster County, Nebraska. We took the mean week end percent harvested of these two crops during the study season and compared these to the rural recorded vibratory noise levels. The week end percent harvested was positively correlated with the daily average Leq for rural sites (Chisq = 4.4975, P = 0.034, conditional R^2 = 0.381).

24 Hours

## [1] 2637 2661
## quartz_off_screen 
##                 2

## quartz_off_screen 
##                 2

## Linear mixed model fit by REML ['lmerMod']
## Formula: meanleq ~ dark_light * Category + (1 | Site)
##    Data: houravgcatl
## 
## REML criterion at convergence: 35380.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0671 -0.6457 -0.1026  0.5184  5.9609 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 8.152    2.855   
##  Residual             8.497    2.915   
## Number of obs: 7080, groups:  Site, 23
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)              -60.239772   0.696933 -86.435
## dark_light                -0.170665   0.005835 -29.248
## CategoryRural             -5.747389   1.364398  -4.212
## dark_light:CategoryRural  -0.004469   0.011348  -0.394
## 
## Correlation of Fixed Effects:
##             (Intr) drk_lg CtgryR
## dark_light  -0.096              
## CategoryRrl -0.511  0.049       
## drk_lght:CR  0.050 -0.514 -0.096
##            R2m       R2c
## [1,] 0.3233593 0.6546749

## quartz_off_screen 
##                 2

Here we assessed how vibratory noise levels change throughout the day. We used hours since sunrise since preliminary analysis revealed that vibratory noise is likely highest at dawn. From dawn, vibratory noise decreases through the rest of the day and night (Chisq = 1179.0765, P < 0.001, conditional R^2 = 0.655). Still, vibratory noise significantly differs by category, with urban sites being louder than rural sites (Chisq = 18.230, P < 0.001).

## quartz_off_screen 
##                 2

## quartz_off_screen 
##                 2

## quartz_off_screen 
##                 2